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ABSTRACT 

A  time-dependent,  two-dimensional  thermodynamic  model  of  an  open  lead 
in  central  arctic  sea  ice  is  presented.  The  model  is  generated  by  opening 
a  lead  of  finite  width  and  infinite  extent  in  the  equilibrium  sea-ice  cover. 
From  this  initial  condition,  the  model  is  integrated  numerically  as  a  sea- 
ice  cover  is  reestablished  over  the  lead.  The  effects  of  various  representa- 
tive advective  parameterizations,  and  temperature  and  salinity  profiles  in 
the  ocean's  surface  layer,  on  the  heat  flux  through  the  lead  and  on  the 
nature  of  ice  formation  are  investigated  as  the  lead  closes  by  thermal 
processes  alone.  Continuity  equations  involving  a  horizontal  advection 
term  and  a  vertical  diffusion  term  govern  heat  and  salt  transport  in  the 
water.  Cooling-induced  convective  overturn  as  a  mechanism  for  vertical 
heat  and  salt  transport  in  the  water  column  is  treated  through  an  artifice 
of  the  vertical  diffusion  term. 
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I.     INTRODUCTION 

A.      POLAR  AIR-SEA  INTERACTIONS 

Heat  lost  from  polar  regions  together  with  the  heat  source  of  the 
tropics  form  the  basis  of  the  global   general   circulation.     The  flux  of  heat 
into  polar  regions  from  lower  latitudes  results  from  the  general   meridional 
temperature  gradient  from  polar  regions  toward  the  equator  due  to  unequal 
latitudinal   heating  of  the  earth's  surface.     Most  of  the  heat  advected  into 
polar  regions  by  the  atmosphere  is  lost  to  space  by  long  wave  radiation, 
thus   forming  a  basic  mechanism  for  maintaining  planetary  radiative 
equilibrium. 

Heat  is  also  advected  into  polar  regions  by  ocean  currents,  but  due 
to  the  relative  inefficiency  of  air-sea  heat  transfer  and  the  relative 
slowness  of  ocean  currents,   this  transfer  is  roughly  one-tenth  that 
attributed  to  the  atmosphere. 

Sea  ice  is  a  dominant  feature  of  both  arctic  and  antarctic  seas;  however, 
due  to  major  differences  in  the  physical   configuration  of  sea  and  land 
masses  of  these  regions,   it  plays  decidedly  different  roles.     The  arctic, 
being  essentially  an  ice-covered  ocean,   is  extremely  sensitive  to  climatic 
perturbations  since  the  areal   extent  and  depth  of  sea  ice  is  critical   to  the 
regional   energy  balance.     A  surplus  of  about  one-third  of  the  annual   energy 
influx  would  remove  the  arctic  ice  pack   [Maykut  and  Untersteiner  1969,  1971  J. 
On  the  other  hand,   the  antarctic  continental   mass,  covered  by  2000  to  4000 
meters  of  ice,   is  relatively  insensitive  to  small   energy  flux  perturbations, 
although  its  sea-ice  cover  would  be  affected  similarly  to  that  of  the  arctic. 


Of  major  import  to  polar  air-sea  interactions  is  the  degree  to  which 
the  sea  surface  is  covered  by  ice.  Sea  ice  acts  analogously  to  a  'lid' 
on  the  sea,  modifying  the  dynamics  and  continuity  of  the  sea  and  the  air- 
sea  system.  Sea-ice  cover  dramatically  reduces  the  heat  transfer  between 
ocean  and  atmosphere  by  effectively  limiting  that  transfer  to  a  process  of 
molecular  conduction.  Further,  by  intercepting  the  vertical  flux  of 
momentum  transferred  from  the  mean  wind  field  to  the  sea  surface  by  the 
Reynold's  stress,  it  suppresses  wind  mixing.  Finally,  it  reflects  a 
considerable  portion  of  the  incoming  short-wave  radiation,  imposes  an 
upper  limit  on  the  sea's  surface-boundary  temperature,  reduces  latent  heat 
transfer  by  impeding  evaporation,  and,  through  its  changes  of  phase,  serves 
as  a  buffer,  or  damping  mechanism.  Alternatively,  heat  transfer  processes 
are  significantly  enhanced  over  open-water  areas  in  proximity  to  large  ice 
masses,  particularly  during  the  winter,  as  surface  air  temperatures  drop 
far  below  surface  water  temperatures.  Wind  and  water  stresses,  in  addition 
to  directly  affecting  ocean-atmosphere  heat  transfer  processes,  keep  the 
ice  in  nearly  continuous  motion,  thus  locally  affecting  heat  transfer  as 
well  as  producing  mechanical  stresses  resulting  in  the  formation  of  pressure 
ridges,  polynya,  and  open  leads. 

Polar  air-sea  interactions  are  therefore  characterized  by  complex  feed 
back  wherein  the  ocean  and  atmosphere  individually  and  collectively  influence 
the  thickness  and  areal  extent  of  the  ice,  which  in  turn  has  a  pronounced 
effect  on  the  exchange  in  question. 

B.  SURVEY  OF  PREVIOUS  RESEARCH 

Range  and  depth  of  interest  in  arctic  research  has  expanded  considerably 
since  Nansen's  epic  drift  in  the  FRAM  (1893-1896).  From  oioneering  studies 


until  only  recently,  research  emphasis  has  been  primarily  of  a  descriptive 
nature,  and  it  was  not  until  the  Soviet  establishment  of  the  drifting 
station  "North  Pole  1  (NP-1)"  in  1937  that  a  systematic  study  of  the 
arctic  interior  was  initiated.  Launching  of  a  comparable  long-range 
arctic  research  program  by  the  United  States  would  await  the  establishment 
of  "Fletcher's  Ice  Island  (T-3)"  in  1952. 

Attempts  to  understand  the  response  of  sea  ice  to  environmental 
changes  has  long  stimulated  the  interest  of  arctic  researchers.  Quantita- 
tive prediction  of  these  responses  has  classically  involved  two  primary 
approaches,  the  empirical  and  the  theoretical.  Although  providing  minimal 
insight  to  the  dynamics  of  physical  processes  occurring  between  air,  ice, 
and  sea,  empirical  formulae  relating  sea-ice  growth  processes  to  observed 
surface  temperatures  have  accomplished  basic  predictive  goals.  The  most 
sophisticated  example  of  this  approach  is  perhaps  one  advanced  by  Bilello 
[1961 j.  Theoretical  studies  involved  analytical  solutions  and  were  there- 
fore severely  limited  in  scope  and  complexity.  One  of  the  foremost  of 
these  analytical  treatments  is  attributed  to  Kolesnikov  [U.S.  Naval  Civil 
Engineering  Laboratory  1966]. 

In  the  early  1950' s  the  tenor  of  arctic  research  began  to  reflect 
increased  emphasis  on  the  theoretical  approach  to  arctic  dynamics  as 
comprehensive  data  from  systematic  studies  became  available.  It  was  at 
this  point  that  studies  relating  the  arctic  energy  balance  to  the  mass 
budget  and  subsequent  research  directed  toward  heat  transfer  processes 
began  to  emerge.  Development  of  physical  relationships  was  confined  to 
empirical  and  statistical  approaches,  or  theoretical  problems  to  which 
simple,  analytical,  approximate  solutions  could  be  developed. 
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Classic  arctic  heat  budget  studies  were  presented  by  several  authors. 
Badgley  [1961,  1966]  summarized  and  collated  available  data  to  deduce  a 
typical  heat  budget  at  the  surface  of  the  central  Arctic  Ocean.  Fletcher 
[1966]  discussed  annual  patterns  of  atmospheric  heat  loss  for  each  component 
of  the  heat  budget  in  relation  to  the  general  atmospheric  circulation,  for 
both  ice-covered  oceans  and  ice-free  oceans.  Vowinckel  [Arctic  Meteorology 
Research  Group  1964]  independently  investigated  the  amount  of  heat  conducted 
through  the  Arctic  Ocean  ice  to  check  previously  calculated  values  for  the 
heat  released  from  the  Arctic  Ocean.  Doronin  [1966]  presented  an  analysis 
of  heat-balance  components  of  the  surface  layer  in  the  Arctic  Ocean. 
Coachman  [1966]  proposed  a  model  to  explain  the  observed  features  of  the 
oceanic  regime  of  the  Arctic  surface  layer  during  the  winter.  Muench 
[Baffin  Bay  -  North  Water  Project  1971]  discussed  the  physical  oceanography 
of  the  northern  Baffin  Bay  region  in  considerable  detail. 

It  was  not  until  the  1960's  and  the  attendant  widened  use  of  the  high- 
speed digital  computer  as  a  research  tool  that  physical  processes  could  be 
modeled  free  of  constraints  imposed  by  analytical  solutions.  Using  finite 
difference  schemes,  a  system  of  integral-differential  equations  could  now 
be  solved  rapidly  within  acceptable  error  limitations.  Typical  of  first- 
generation  numerical  models  in  the  field  is  that  developed  by  Untersteiner 
[1966],  which  predicted  thickness  and  temperature  of  sea  ice,  thus  allowing 
further  insight  into  the  relationship  of  climatic  change  to  ice  production, 
ice  decay,  and  equilibrium  ice  thickness.  The  advanced  one-dimensional 
thermodynamic  model  of  central  Arctic  sea  ice  developed  by  Maykut  and 
Untersteiner  [1969,  1971]  is  a  prime  example  of  the  present  theoretical 
approach. 
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C.  PRESENT  RESEARCH  DIRECTION 

A  great  deal  of  research  effort  is  presently  being  directed  toward 
developing  numerical  models  of  the  atmosphere  and  of  the  ocean  as  well 
as  atmosphere-ocean  models,  with  the  objective  of  ultimately  producing 
an  efficient,  reliable  global  system  model.  Problems  inherent  in  the 
realization  of  this  objective  are  numerous  and  include  not  only  a  lack  of 
refinement  in  respective  models  of  the  atmosphere  and  ocean,  but  also  a 
lack  of  refinement  in  the  technique  of  coupling  the  respective  models. 
Further,  observational  input  data  is  in  general  deficient  for  these  models, 
particularly  in  the  Arctic.  Present  state-of-the-art  examples  of  models 
of  this  type  are  well  represented  by  the  sophisticated  models  of  Manabe 
[1969]  and  Bryan  [1969]. 

One  of  the  more  ambitious  research  projects  in  recent  years  is  the 
Arctic  Ice  Dynamics  Joint  Experiment.  AIDJEX  is  essentially  a  United 
States  -  Canadian  cooperative  effort  to  gain  quantitative  understanding 
of  the  interaction  between  atmospheric  sea  ice  and  fluid  ocean  fields  of 
motion.  In  progress  at  the  time  of  this  writing,  AIDJEX  is  unique  in  its 
observational  approach  since,  in  order  to  gain  a  full  understanding  of 
interacting  fields,  it  involves  the  collaboration  of  many  research  groups 
in  making  coordinated  measurements  of  interacting  fields  of  motion,  stress, 
and  strain  -  observing  appropriate  time  and  space  scales  -  over  a  minimum 
period  of  one  year.  The  fruits  of  this  ambitious  and  imaginative  study, 
in  terms  of  fresh  insights,  controlled  data,  and  international  cooperation 
appear  to  be  vast  already. 
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D.  CONCEPTUALLY-RELATED  STUDIES 

Studies  of  Arctic  dynamics,  energy  balances  and  mass  budgets  range 
over  a  wide  spectrum.  Scientific  literature  reflects  numerous  experimental 
and  theoretical  models  of  microstructure  and  micro-scale  energy  transfer 
processes  between  ocean  and  atmosphere,  ice  and  atmosphere,  and  ice  and 
ocean,  essentially  occurring  in  boundary-layer  regions.  Pertinent  to  the 
work  presented  herein  is  research  performed  by  Foster  [1968,  1969]  in  which 
haline  convection  induced  by  the  freezing  of  sea  water  was  studied,  and  by 
Lake  and  Lewis  [1970]  who  studied  salt  rejection  by  sea  ice  during  growth. 
Work  performed  by  Lewis  and  Walker  [1970]  illustrates  seasonal  changes  in 
the  temperature  and  salinity  profiles  beneath  annual  sea  ice  cover,  and 
Stern  and  Turner  [1969]  reported  on  the  formation  of  salt  fingers  and 
convecting  layers. 

Directly  related  to  the  model  proposed  herein  is  the  macroscale  advanced 
model  derived  by  Maykut  and  Untersteiner  [1969,  1971]  to  forecast  ice 
temperature  and  thickness,  which  defines  the  role  of  each  component  of  the 
energy  budget  in  relation  to  its  effects  upon  the  ice.  Although  comprehen- 
sive in  its  treatment  of  physical  processes  related  to  components  of  the 
energy  budget,  the  model  admittedly  has  several  shortcomings.  Physical 
limitations  are  due  to  uncertainties  in  environmental  data  and  include 
the  necessary  assumption  of  a  constant  oceanic  heat  flux  and  generally 
course  albedo  and  energy-flux  data.  Limitations,  arising  from  an  inability 
to  model  certain  processes  include:  1)  the  neglect  of  mechanical  stresses 
produced  by  winds  and  surface  currents,  2)  treatment  of  turbulent  fluxes 
at  boundaries  as  being  independent  of  growth  rates  or  the  physical  states 
of  ice  and  snow,  3)  imprecise  treatment  of  downward  heat  transport  from 
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melting  snow,  4)  undefined  shapes  of  salinity  profiles  for  various  ice 
conditions,  and  5)  neglect  of  heat  storage  by  melt  ponds.  Leads  were 
considered  only  so  far  as  they  indirectly  affected  the  energy  fluxes  over 
the  bare  ice.  In  spite  of  its  limitations,  this  will  certainly  be  a 
cornerstone  for  future  thermodynamic  ice  models  and  provides  a  number  of 
conceptual  bridges  used  in  the  present  work. 

E.  THESIS  OBJECTIVE 

It  is  known  that  significant  heat  transfer  occurs  between  the  atmosphere 
and  the  ocean  through  open  leads  in  the  Arctic  sea  ice  cover.  As  found 
experimentally  by  Badgley  [1961],  sensible  heat  loss  is  at  least  two  orders 
of  magnitude  greater  from  open  leads  than  from  perennial  sea  ice. 

With  the  need  for  further  refinement  in  models  of  the  ocean  and  the 
atmosphere,  as  well  as  in  model  coupling  schemes,  the  requirement  for 
mesoscale  examination  of  ocean-ice-atmosphere  heat- transfer  processes  is 
immediate  and  apparent.  With  appropriate  observational  input  data  becoming 
increasingly  more  abundant,  informative  numerical  models  of  mesoscale 
processes  appear  promising  for  furthering  the  understanding  of  physical 
processes  and  for  realizing  control  of  mass-energy  fluxes  in  the  arctic. 

The  objective  of  this  study  was  to  develop  a  thermodynamic  model  of  an 
open  lead  in  central  arctic  sea-ice  in  which  hydrodynamic  processes  of 
convection,  advection,  and  diffusion  of  heat  and  salt  were  represented. 
By  integrating  the  model  numerically,  the  effects  of  various  representative 
advective  parametizations,  and  temperature  and  salinity  profiles  in  the 
ocean's  surface  layer,  on  the  heat  flux  through  the  lead  and  on  the  nature 
of  ice  formation,  could  be  investigated  as  the  lead  reclosed  by  thermal 
processes  alone. 
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II.   NATURE  OF  THE  PROBLEM 

Since  sea  ice  essentially  insulates  the  ocean  from  the  atmosphere  and 
suppresses  important  thermodynamic  processes,  it,  along  with  heat-advecting 
ocean  currents,  becomes  a  primary  factor  in  determining  the  vertical 
structure  in  Arctic  seas.  However,  the  dynamic  transition  from  ice-free  to 
ice-covered  seas  is  not  fully  understood.   In  particular,  the  nature  of 
the  fluxes  of  salt  and  heat,  as  well  as  the  structure  and  depth  of  unstable 
convective  overturn  during  initial  phases  of  the  freezing  process,  require 
study. 

A.  DYNAMICS  OF  THE  SEA  WATER  FREEZING  PROCESS 

The  structure  of  the  active  surface  layer  of  the  Arctic  Ocean  is 
generally  horizontally  uniform  with  a  very   stable  density  gradient.  How- 
ever, during  the  freezing  process  the  density  gradient  may  be  significantly 
altered  by  three  primary  mechanism:  1)  progressive  cooling  of  surface 
waters  from  above,  2)  brine  expulsion  by  the  ice  during  the  freezing  pro- 
cess, and  3)  advection  of  heat  and/or  salt  by  surface  currents. 

The  fundamental  physical  relation  which  makes  the  sea  water  freezing 
process  unique  may  be  understood  by  examining  the  curves  of  maximum  density 
and  freezing  point  for  sea  water,  plotted  with  respect  to  temperature  versus 
chlorinity  or  salinity  (Figure  1).  Specifically,  for  waters  of  salinity 
greater  than  24.70  parts  per  thousand,  the  temperature  of  maximum  density  of 
the  liquid  state  will  always  be  lower  than  the  liquid's  freezing-point 
temperature  under  pressures  existing  in  nature.  Any  process  that  tends  to 
produce  a  vertical  density  instability  will  induce  convective  overturn  in 
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the  water  column,  resulting  in  mixing  of  the  unstable  segment  such  as 
would  be  produced  in  the  freezing  process  by  surface  cooling  and  freezing- 
associated  brine  expulsion.  Convective  overturn  will  therefore  be  continuous 
in  the  unstable  segment  of  the  liquid  water  column  from  the  onset  of  surface 
cooling  to  the  point  of  total  change-of -phase  of  the  water  column.  A  some- 
what less  obvious  but  straightforward  point  should  be  mentioned  in  qualifi- 
cation of  the  foregoing.  From  an  arbitrary  water  temperature  above  the 
freezing  point,  down  to,  but  not  including  the  freezing  point,  convection 
is  temperature  controlled  during  cooling  since  there  is  no  ice  formed, 
hence,  no  brine  expulsion.  During  cooling  of  surface  waters  at  the  freezing 
point  for  existing  surface  salinities,  convection  is  primarily  salinity 
controlled  since  brine  is  being  continuously  expelled  from  newly  forming 
ice.  Within  this  latter  water  temperature  regime,  temperature  changes  play 
a  minor  role  in  effecting  changes  in  surface  water  densities. 

As  convective  overturn  is  a  mechanism  for  inhibiting  freezing  by 
continuously  transferring  heat  from  depth  to  the  active  boundary  at  which 
freezing  is  taking  place,  so  it  is  also  a  mechanism  for  effectively  removing 
heat  from  depths  significantly  greater  than  could  be  effected  by  molecular 
diffusion.  An  equilibrium  ice  thickness  is  reached  when  the  resultant  heat 
flux  upward  through  convection,  advection  and  diffusion  in  the  water  column 
balances  the  upward  resultant  heat  flux  through  the  ice,  at  the  ice-water 
boundary.  Prior  to  this  point,  ice  is  accreted  in  depth  since  the  upward 
heat  flux  through  the  ice  exceeds  the  upward  heat  flux  in  the  water  column 
arriving  at  the  ice-water  boundary. 
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B.  DISTRIBUTION  OF  VARIABLES 

The  distribution  of  temperature,  salinity  and  velocity  are  considered 
within  a  volume  of  water  under  a  refreezing  lead.  The  water  is  assumed  to 
be  inviscid,  horizontally  homogeneous,  horizontally  isotropic,  irrotational , 
and  incompressible.  Accelerations  resulting  from  the  coriolis  effect  and 
horizontal  pressure  gradients  are  assumed  to  be  negligible.  Momentum,  heat 
energy,  and  salt  mass  are  conserved  within  the  volume  such  that  no  velocity, 
temperature,  or  salinity  sources  or  sinks  exist  except  at  the  volume's 
boundaries.  Furthermore,  currents  are  assumed  to  be  steady-state  and 
uniform  both  horizontally  and  vertically. 

Using  the  classic  box-model  approach,  these  conservative  fields  are 
represented  continuously  within  the  volume  by  individual  continuity 
equations.  For  the  vector  field  of  velocity,  where  K  is  the  diffusivity 
and  Sw  is  a  source  term. 

DV/_  =  V-KW+  S  (1) 

Dt  v 

which  reduces  to 

v(t,x,y,z)  =  v  =  constant  (2) 

i.e.,  the  velocity  is  specified  as  uniform  and  steady. 

Similarly,  the  scalar  fields  of  heat  and  salt,  here  collectively 
depicted  by  the  scalar  quantity  0,  are  represented  by 

DJ3/Dt  =    v.    K  v  0  +  S  (3) 


whe 


re  K  is  the  diffusivity  and  S.  is  a  source  term. 
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This  reduces  to 

*%t  =  "  V'V(3  +  Khv20  +Vv   '   K(z)  vv0  (4) 

where  K  is  the  horizontal  diffusivity  and  vv  is  the  vertical  del  operator 
since, 

K(x)  =  K(y)  =  KH  =  constant 

The  local,  or  time-rate-of-change,  is  therefore  equal   to  the  sum  of  the  field- 
rate-of-change  and  the  diffusive-rate-of-change. 

C.     HYDRODYNAMIC  PROCESSES 

Mechanisms  for  the  hydrodynamic  transport  of  heat  and  salt  within  the 
lead  are  advection  and  diffusion.     However,  from  the  dynamics  of  the  sea- 
water  freezing  process,   it  is  known  that  convection  plays  a  major  role  in 
transporting  both  heat  and  salt  in  the  water  column.     Although   it  is  riot 
represented  by  a  discrete  term  in  the  foregoing  continuity  equation,  the 
convective  transport  term  may  be  thought  of  as  simply  a  temporary  (when 
vertical   density  instabilities  arise)  modification  of  the  vertical   diffusion 
term  taking  the  form 


'im  *A       «■)""'*,  (5) 


K(*)-« 

which  in  essence  results  in  free  vertical  (convective)  transport. 

D.  MASS  AND  ENERGY  BALANCE 

With  the  foregoing  analyses  of  continuity  relations  and  mechanisms  of 
hydrodynamic  transport  existing  within  the  box  model  of  the  open  lead,  the 
mass  and  energy  balances  existing  at  the  volume's  boundaries  may  be 
examined. 


19 


Over  an  arbitrary  surface  area  of  a  reclosing  lead,  a  given  change  of 
heat,  AQ,  at  the  water  surface  must  equal  the  heat  added  by  advection  and 
diffusion  at  upstream  boundaries,  Q  ,  plus  the  heat  added  from  depth  in  the 
lead  by  convection,  advection,  and  diffusion,  Qz,  plus  the  heat  added  by 
surface  ice  accretion,  Qf,  plus  the  open-water  heat  loss  at  the  surface  by 
conduction,  evaporation  and  radiation,  Q  +Q  +Q  ,  plus  the  heat  loss  through 
the  ice  cover  by  conduction,  Q. ,  minus  the  heat  diffused  across  vertical 
boundaries  normal  to  advective  transport,  Qn,  minus  the  heat  transported 
across  the  downstream  boundaries  by  advection  and  diffusion,  Q^,  such  that 
AQ  =  Qu  +  Qz  +  Qf  +  W(QS+Qe+Qr)  +  Qi (1-W)  -  Qn  -  Qd  (6) 

Where  W  is  the  fraction  of  open  water  comprising  the  surface  of  the  lead  at 
any  given  time.  The  heat  added  from  depth  in  the  lead  is  heat  ultimately 
derived  from  the  Bering  Sea  Water  (in  the  Canadian  Basin)  and  from  the 
Atlantic  layer. 

By  similar  analysis,  a  given  change  of  salinity,  aS,  at  the  water's 
surface  must  equal  the  salt  added  by  advection  and  diffusion  at  upstream 
boundaries,  Su,  plus  the  salt  added  from  depth  in  the  lead  by  convection, 
advection,  and  diffusion,  Sz,  plus  the  salt  added  by  surface  ice  accretion 
(salt  rejection),  Sf,  plus  the  salt  lost  through  surface  ice  melting  (salt 
dilution),  Sm,  plus  the  salt  added  by  surface  evaporation,  Se,  minus  the 
salt  diffused  across  vertical  boundaries  normal  to  advective  transport, 
Sn,  minus  the  salt  transported  across  downstream  boundaries  by  advection  and 
diffusion,  S^,  such  that 

AS  =  su  +  sz  +  (i-w)(sf+sm)  +  w(se)  -  sn  -  sd  (7) 

where  again  W  is  the  fraction  of  open  water  comprising  the  surface  of  the 
lead  at  any  given  time.  The  salt  added  from  depth  in  the  lead  is  derived 
from  the  generally  strong  positive  halocline.  existing  at  depths  below  25  to 
50  meters. 
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Qualitatively  then,  before  the  equilibrium  sea  ice  cover  has  been 
reestablished,  there  will  exist  differentials  in  the  temperature  and 
salinity  profiles  of  water  columns.  A  primary  process  of  interest  in  the 
refreezing  lead  is  that  of  surface  ice  formation,  which  lends  a  unique 
character  to  lead  hydrodynamics.  The  heat  added/lost  by  surface  ice 
accretion/melting  is  a  function  of  the  incremental  change  in  the  mass  of 
ice,AI,  and  the  heat  of  formation  of  sea  ice,  L . ,  such  that 

Qf  =  -Qm  =  L-AI  (8) 

Coupled  directly  to  this  process  is  the  process  of  brine  expulsion  (salt 
rejection)  by  the  sea  ice  during  growth.  The  salt  added/lost  by  surface  ice 
accretion/melting  is  a  function  of  the  incremental  change  in  the  mass  of  ice, 
AI,  the  salinity  of  the  sea  ice,  S^ ,  (which  in  turn  is  dependent  upon  rate- 
of-freezing,  ice  temperature,  and  age  of  the  ice),  and  the  salinity  of  the 
parent  sea  water,  Sw,  such  that 

Sf  =  -Sm  =  f(Al,Si,Sw)  (9) 

Qualitative  analysis  of  the  salt  rejection  process  by  constituent  salts  is 
not  appropos  "to  either  the  scale  or  the  nature  of  the  study  at  hand;  there- 
fore salt  rejection  was  treated  simply  as  a  process  of  water-salinity 
modification. 
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III.   EXPERIMENTAL  PROCEDURE 

A.  FORMULATION  OF  THE  MODEL 

1 .  Establishment  of  a  Representative  Cross  Section 

In  deciding  upon  what  physical  dimensions  would  be  appropriate  for 
the  model,  it  appeared  most  suitable  to  represent  a  small  open  lead  in 
central  arctic  sea  ice  such  as  one  would  expect  to  observe  between  early 
to  late  winter.  Therefore,  the  width  of  the  lead  was  chosen  to  be  150 
meters,  with  an  additional  50  meters  under  equilibrium  sea-ice  cover  on 
the  downstream  side  of  the  model  lead  to  facilitate  studying  resultant 
thermodynamics  and  hydrodynamics  in  that  region. 

The  model's  depth  extent  was  set  at  50  meters  based  upon  lower 
boundary  condition  considerations.  Over  much  of  the  Arctic  Basin  the 
vertical  temperature  gradient  at  that  depth  is  zero  or  yery   slightly  posi- 
tive. Since  a  no-flux  lower  boundary  condition  was  desirable  in  order  to 
better  study  the  effects  of  advection  of  conservative  constituents  in  the 
lead,  that  depth  was  chosen. 

It  was  felt  that  the  model  could  be  simplified  without  the  loss  of 
any  essential  physics  if  derivatives  along  the  lead  (Y-axis)  were  set  to 
zero.  This  reduced  the  problem  to  one  of  looking  at  a  cross-sectional 
distribution  in  the  lead.  Figure  2  is  a  pictorial  depiction  of  the  open 
lead. 

2.  Adaptation  of  the  Continuity  Equations  to  the  Model 

Using  the  simplification  of  the  previous  paragraph  which  allows  the 
use  of  a  cross-sectional  area  rather  than  a  true  volume  to  represent  the 
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Figure  2. 

Pictorial   Depiction  of  Open  Lead 
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lead,  the  horizontal  advection  and  horizontal  diffusion  terms  describing 
changes  along  the  Y-axis  were  eliminated  from  equation  (4). 

The  problem  was  further  simplified  by  assuming  that  all  horizontal 
diffusion  was  negligible  relative  to  horizontal  advection,  and  that  verti- 
cal advection  could  be  disregarded. 

Applying  the  foregoing  simplifications  to  the  expanded  form  of 
equation  (4),  and  choosing  the  state  variables  of  temperature  and  salinity 
as  the  relevant  dependent  variables,  the  final  form  of  the  equations  used 
to  describe  the  conservative  scalar  fields  of  temperature  and  salinity 
within  the  lead  was 

d0/dt  =  -  udp/ax  +  a/az   K(*)d0/Sz  do) 

where, 

0  =  0  (t,x,z) 

3.  Initial  Conditions 

Sea-water  temperatures  and  salinities  within  the  lead  at  the 
instant  that'the  lead  is  generated  are  represented  by  local  stable  tem- 
perature and  salinity  profiles  previously  existing  beneatii  the  equilibrium 
sea-ice  cover.  The  temperature  and  salinity  fields  are  both  horizontally 
uniform  and  vary  only  with  depth  in  the  lead  such  that 

T(0,X,Z)  =  T(Z)  (11) 

S(0,X,Z)  =  S(Z)  (12) 

The  surface  of  the  newly-generated  lead  is  initially  taken  to  be 
ice-free  over  the  left  (upstream)  150  meters,  and  ice-covered  with  equili- 
grium  thickness  ice  over  the  right  (down-stream)  50  meters.  Surface 
temperatures  are  initially  at  the  freezing  point,  in  temperature  equilibrium 
with  the  preexisting  sea-ice  cover. 
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4.  Boundary  Conditions 

The  model's  left,  or  upstream  boundary,  may  be  regarded  as  a 
source  boundary  across  which  sea  water  having  a  given  temperature  and 
salinity  profile  flows  into  the  lead.     Sea-water  temperature  and  salinity 
profiles  are  taken  as  steady-state  functions  of  depth  alone,  identical   to 
the  lead's   initial   temperature  and  salinity  profiles  such  that 

T(t,0,Z)  =  T(Z)  (13) 

S(t,0,Z)  =  S(Z)  (14) 

The  model's  surface  boundary  is  the  primary  region  of  mass  and 
energy  fluxes  and  is  therefore  the  most  dynamically  active  boundary.     The 
boundary  condition  on  temperature  across  the  entire  surface  is  that  of  an 
upward  (or  downward)  diffusive  flux  of  heat  across  the  water  surface  such 
that 

%x  T(t,X,0)       =  f(AT,Zi,Si,Kw,Ki)  (15) 

where  aT  is  the  air-water  temperature  difference   (an  externally  specified 
parameter),  Z-j   is  the  surface  ice  thickness,  S-   is  the  surface  ice  salinity 
(an  externally  specified  constant  parameter),  and  Kw  and  K-j   are  the  co- 
efficients of  thermal   eddy  diffusivity  of  sea  water  and  thermal   conducti- 
vity of  sea  ice  respectively.     When  no  ice  is  present  at  the  surface,  the 
surface  heat  flux  is  taken  as  the  cumulative  open-water  sensible,  latent, 
and  radiative  heat  fluxes,  which  are  specified  in  the  model   as  external 
parameters  according  to  the  air-water  temperature  difference.     Since  the 
model's  surface  at  all   times  represents  the  water  surface,  the  fluxes  at 
that  boundary  are  matched  as  ice  is  formed,  such  that  the  flux  through  the 
water  surface  is  always  equal   to  the  flux  through  the  ice  to  the  atmosphere, 
but  at  no  time  does   it  exceed  the  maximum  open-water  flux  to  the  atmosphere. 
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The  boundary  condition  on  sea  water  salinity  at  the  surface 
specifies  a  no-diffusive  flux  condition  between  air  and  water  as  well  as 
between  ice  and  water  and  ice  and  air  such  that 

%t     S(t.x.o)   =  0  (16) 

However,  salt  rejection  from  ice  to  water  during  the  freezing  process   (as 
well   as  salt  dilution  during  melting)   is  represented  as  a  discrete  process 
of  physically  adding  the  change  in  salinity  produced  by  sea-ice  accretion 
or  melting  to  the  surface  waters  such  that 

S(t,X,o)=»    S(t,X,o)    +  f(Al,AS)  (17) 

where  AI  is  the  incremental  change  in  the  mass  of  sea  ice,  and  AS  is 
the  salinity  differential  existing  between  the  sea  ice  and  its  parent  sea 
water. 

The  bottom  boundary,  as  previously  discussed,  is  characterized 
as  a  no-flux  boundary,  thus 

*/^m   T(t,X,h)    =  0  (18) 

%t    S(t,X,h)         =  0  (19) 

where  h  is  the  model's  depth  extent. 

The  right,  or  downstream  boundary,   is  also  characterized  as  a 
no-diffusive  flux  boundary  where  whatever  temperature-salinity  field  is 
generated  by  the  thermodynamic  processes  of  the  lead  is  simply  advected 
out.     However,  as  mentioned  previously,  the  right  50  meters  of  the  lead 
model's  surface  is  overlain  by  an  equilibrium  sea-ice  cover.     This  boundary 
condition  is  therefore  representative  of  a  return  towards  the  equilibrium 
vertical   profiles.     At  the  model's  right  boundary  then 
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a/dX     T(t,W,Z)       =  0  (20) 

^d  X     S(t,W,Z)       =  0  (21) 

where  W  is  the  model's  width. 

Figure  3  describes  processes  considered  in  the  model   cross  section. 
Figure  4  schematically  represents  the  model   cross  section  with  its  governing 
continuity  equations  and  boundary  conditions.      In  this  figure,  the  surface 
boundary  condition  is  simply  a  statement  of  the  balance  of  fluxes  at  the 
ice-water  boundary.     With  no  ice  cover  (5=0),  the  upward  heat  flux  through 
the  water  surface,   Fz    is  equal   to  the  sum  of  open-water  sensible   (F  ),   latent 
(F  ),  and  radiative   (F  )  heat  fluxes.     When  an  ice  cover  is  present   (5>0),  the 
upward  heat  flux  through  the  water  surface   (*  =0+)   is  equal   to  the  upward  heat 
flux  through  the  ice   ( *  =0")  minus  the  heat  added  by  surface  ice  accretion. 
The  above  relations  are  presented  in  terms  of  temperature  gradients  in  Figure 
4  as  they  are  applied  to  the  governing  continuity  equation  for  temperature. 
Symbols  not  previously  defined  are  p.,  volume  density  of  sea  ice;  q . ,  heat  of 
formation  of  sea  ice;    5  ice  thickness   (volume  of  ice  per  unit  surface  area). 
5.     The  Finite-Difference  Scheme 

A  purely  explicit  finite  difference  scheme  was  selected  to  numeri- 
cally represent  the  linear  parabolic  partial   differential   equation   (10). 
The  time  rate-of -change  term  was  represented  by  forward  differencing,  the 
field  rate-of-change  term  by  backward   (upstream)   differencing,  and  the 
diffusive  rate-of  change  term  by  central   differencing.     Thus  over  the  interior 
(non-boundary)   region  of  the  model 

Ck  -  "j.r  . .  »j.k  -»j-i,k  +  K  g",k+r2gj/gj,k-i    m 

At  AX  (AZ)2 
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Figure  3. 

Description  of  Processes  Considered  in  Model 
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(1)  Advection  of  steady-state  temperature-salinity  field  into  the 
lead  from  under  equilibrium  ice  pack. 

(2)  Heat  loss  to  the  atmosphere  as  a  function  of  ice  thickness, 
ice  accretion/melting,  and  salt  rejection/dilution. 

(3)  Heat  and  salt  transport  by  convective  overturn,   advection  and 
eddy  diffusion. 

(4)  No  transport  across  boundary. 

(5)  Advection  of  temperature-salinity  field  generated  by  lead 
thermodynamics  out  of  the  lead  model. 
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where, 

0n    =  0nAt  ,  (23) 

■j,k     JAX,  kAZ  ^6> 

At  the  left  boundary,  equation  (22)  was  inapplicable  since  tempera- 
ture and  salinity  were  both  steady-state  functions  of  depth  alone.  There- 
fore, at  this  boundary  at  all  times 

C  -  0S,k  (24) 

At  the  surface  boundary,  the  diffusive  term  for  temperature  was 
modified  in  equation  (22)  since 


Tn  -Tn 

'j,k+l      'j,k-l 

-    =   ?    T  25 


2AZ 


or, 


since 


•j,k-i  -  'j,k+i     2AI^v  n  U&J 

The  diffusive  salt  term  at  the  surface  was  modified  in  equation  (22) 


Sn  -  sn 

J>k+1         J'""1     =  0  (27) 

2  A  I 


or, 


Sj  ,k-l   =  sj,k+l  (28) 

At  the  bottom  boundary,  the  diffusive  rate-of-change  terms  of  both 
temperature  and  salinity  were  modified  in  equation   (22)   since 

n  an 

Pj,k+1    "  "j,k-l        _  (29) 


2AZ 


=    0 
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or 


0j,k+l  -  »j.k-l 


(30) 


Finally,  at  the  right  boundary,  equation  (22)  remains  unmodified 
since  advection  is  the  only  transport  mechanism  operating  across  it. 

Stability  analysis  of  the  finite  difference  equation  using  the 
Fourier  series  method  [Smith  1965]  shows  that 


n+1 


s.k[»B-2*y 


(31) 


where  e   represents  an  arbitrary  error  introduced  into  the  problem.     There- 
fore, the  finite-difference  equation   (22)   is  stable  when  the  absolute  value 
of  the  amplification  factor  remains  less  than  or  equal   to  unity 


Ati 


U 


K 


1 


<  1 
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6.  Representation  of  Physical  Processes 

Convective  overturn  is  treated  nominally  through  an  artifice  of  the 
vertical  diffusion  term  as  in  equation  (5).  Numerically  however,  convective 
overturn  and  its  associated  mixing  is  accomplished  through  arithmetic 
averaging  of  temperatures  and  salinities  within  regions  of  vertical  density 
instability  as  determined  by  sigma-t  values,  followed  by  recomputation  of 
the  resultant  sigma-t  value  for  the  newly  mixed  (and  neutrally  stable) 
region.  Subroutine  OVRTRN  accomplishes  convective  overturn  in  the  model. 

Advection  of  heat  and  salt  (temperature  and  salinity)  in  the  lead  is 
driven  by  a  current,  u,  in  the  positive-x  direction.  Values  and  depth  pro- 
files of  the  current  as  used  in  the  model  are  given  in  subparagraph  B  of  this 
section. 
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Diffusion  of  heat  and  salt  in  the  vertical  are  determined  by  the 

coefficients  of  vertical  temperature  and  salinity  diffusivities,  both  set 

2     -1 

at  constant  values  of  10  cm  sec    in  the  model. 

Formation  of  sea  ice  is  accomplished  as  in  equation  (8).  More 
specifically  however,  when  a  given  mass  of  water, Aw,  is  cooled  AT  degrees 
centigrade  below  its  salinity  dependent  freezing-point  temperature,  the 
incremental  mass  of  ice  formed,  A  I,  is 

A  I  =  Cp(AW.AT)/,  (33) 

where  Cp  is  the  specific  heat  of  sea  water  and  L-  is  the  heat  of  formation 
of  sea  ice. 

Salt  rejection  (or  dilution)  resulting  from  surface  sea  ice  accretion 

(or  melting)  is  described  in  equation  (17).  Specifically,  if  A I  is  the 

incremental  change  in  the  mass  of  sea  ice  and  AS  is  the  salinity  differential 

existing  between  a  unit  mass  of  parent  sea  water  and  a  unit  mass  of  its 

generated  sea  ice,  then  the  incremental  salinity  change, ASr,  in  the  surface 
waters  resulting  from  that  freezing  or  melting  process  is 

ASR  =  A  I. AS  (34) 

7.     Physical   Quantities  and  Relations 

Sea  water  chlorinity  is  related  to  salinity  [Knudsen  1902]  by  the 

relation 

Cl=(s   -  0.030)/l .8050  (35) 

Freezing-point  depression  is  empirically  related  to  sea  water 

chlorinity  by  Thompson   [Sverdrup,  Johnson  and  Fleming  1942]  by  the  relation 

Tr     =  -0.0966'Cl    -  0.0000052-C13  (36) 

fP 


32 


The  density  of  sea  water  as  a  function  of  salinity  and  temperature 
(sigma-t)   [Knudsen  1902]  is  empirically  computed  through 

*0  =     ^(Cl)  (37) 

AT  =      '2(T)  (38) 

BT  =      '3(T)  (39) 

2T  =     f4(T)  (40) 

't  =     f5(2T»  VAT'BT>  <41> 

The  specific  heat  at  constant  pressure  of  sea  water  is  empirically 
related  to  its  salinity  at  normal    temperature  and  pressure  by  Kuwahara 
[Sverdrup,  Johnson  and  Fleming  1942]  by  the  relation 

C     =   1.005  -  0.004136-S  +  0.0001098-S2 

-  0.000001324-S3  (42) 

The  thermal   conductivity  of  sea  ice  was  determined  using  the  relation- 
ship established  by  Untersteiner  [Maykut  and  Untersteiner   !969,   1971]. 


where, 


ki   =  kf  h-    ^f^     ,  T<0  (43) 

the  subscript  f  refers  to  pure  ice, 

kf  =  0.00486  cal/cm-sec-°C, 

p 
b  =  0.28  cal-cm  /g-sec, 

T  =  temperature  of  ice  in  degrees  Celsius, 

S(z)  =  ice  salinity  at  depth  z  in  g/cm  . 

The  heat  of  formation  of  sea  ice  was  represented  by  the  simplified 

relationship  [Pounder  1965] 


where 


L^Lf.O-f-)  (44) 


Lf  =  79.77  cal/g, 

<*    =  salinity  of  sea  ice,  g/kg, 

S  =  salinity  of  parent  sea  water,  g/kg. 
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B.   ENVIRONMENTAL  DATA  INPUTS 

Environmental  data  used  in  the  model  fell  into  two  basic  categories: 
that  which  remained  constant  throughout  all  runs  of  the  model,  and  that 
which  was  variable  from  run-to-run.  Representative  of  the  former  were 
initial  equilibrium-ice  thickness,  set  at  200  cm,  an  ice  salinity  of  8.0 
g/kg,  and  an  ice  density  of  0.91  g/cm3.  Further,  representative  open  water 
sensible,  latent,  and  radiative  heat  flux  densities  as  presented  by  Muench 
[Baffin  Bay-North  Water  Project  1971]  were  applied  in  the  model.  However, 
as  ice  formed  at  the  surface,  the  latent  heat  flux  density  was  modified  so 
as  to  decrease  exponentially  with  increasing  ice  thickness.  The  radiative 
heat  flux  density  was  modified  so  as  to  decrease  from  open  water  values  as 
ice  formed,  approaching  a  constant  value  representative  of  the  radiative 
flux  density  over  pack  ice  as  also  presented  by  Muench.  The  sensible  heat 
flux  density  was  modified  internally  by  the  relation  shown  in  Figure  4  as 
surface  ice  thicknesses  increased.  All  surface  heat  flux  densities  were 
representative  of  an  air-water  temperature  difference  of  25°C  with  no  wind. 
Finally,  the  coefficients  of  vertical  temperature  and  salinity  eddy  diffusi- 

9 

vity  were  set  at  10.0  cm  /  sec. 

Primary  variable  parameters  in  the  model  were  horizontal  current  magni- 
tudes, vertical  temperature  profiles,  and  vertical  salinity  profiles. 
Currents  applied  in  the  model  were  within  the  range  of  zero  to  10.0  cm/sec. 
Extreme  test  case  values  were  0.0  and  10.0  cm/sec.  Additional  test  case 
values  of  2.0  and  4.0  cm/sec  were  used  to  study  one  particular  temperature- 
salinity  field  combination  with  its  associated  hydrodynamics  in  further 
detail.  A  value  of  7.0  cm/sec  was  chosen  to  represent  a  typical  horizontal 
current  occurring  in  the  lead,  and  was  subsequently  used  as  the  real  case 
value. 
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Test  case  values  of  vertical  temperature  and  salinity  profiles 
representative  of  extremes  in  vertical  gradient  were  applied.  Profiles 
A  represented  strong  gradients,  while  profiles  B  represented  weak 
gradients  (Table  I).  Real  case  temperature  and  salinity  profiles, 
qualitatively  representative  of  early  winter  (profile  C)  and  late 
winter  (profile  D)  conditions  in  the  central  arctic  [Coachman  and 
Barnes  1961,  1962],  were  chosen  to  represent  the  typical  lead. 

Table  I  summarizes  the  profile  designations  of  the  three  variable 
parameters  and  their  various  values. 

C.  MODEL  INTEGRATION  PLAN 

In  order  to  gain  maximum  insight  to  physical  processes  occurring  in 
the  lead,  as  well  as  to  study  the  effects  of  current  magnitude,  vertical 
temperature  gradient,  and  vertical  salinity  gradient  on  lead  thermodynamics, 
trie  i>erieb  01  cesx  Vdi  lauki  was  run  in  nic  mouci  •  rvcprcacHteu  uj   iuh;> 
one  through  eight,  these  test  cases  involved  running  all  possible  combina- 
tions of  current  magnitude  (u)  profiles  A  and  B,  temperature  (T)  profiles 
A  and  B,  and  salinity  (S)  profiles  A  and  B.  These  test  cases  were 
integrated  over  a  lead  reclosure  period  of  two  days. 

As  a  result  of  the  interesting  dynamics  occurring  in  run  6,  additional 
runs  9,  10,  and  14  were  scheduled  using  temperature  and  salinity  profiles 
identical  to  run  6  (and  2),  but  with  current  magnitude  profiles  C,  D,  and 
F  respectively,  in  order  to  further  investigate  these  phenomena.  Integra- 
tion was  again  for  a  period  of  two  days. 

Finally,  the  real-case  variables,  representative  of  naturally 
occurring  parameters,  were  put  into  the  model  and  integrated  over  two-day 
periods  in  runs  11  and  12,  and  over  a  30  day  period  in  run  13. 
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Table  I.     Environmental   Variable  Profiles 


Horizontal 

Current 

Magnitud 

es   (cm/sec) 

Depth 
(m) 

Profile 

A 

B 

C 

D 

E 

F 

0  -  50 

10.0 

1.0 

2.0 

4.0 

7.0 

0.0 

Water  Temperatures  (QC) 


Depth 
(m) 


Profile 


A 


0 

-1.50 

5 

-1.40 

10 

-1.30 

15 

-1.20 

20 

-1.10 

25 

-1.00 

30 

-0.90 

35 

-0.80 

40 

-0.70 

45 

-0.60 

50 

-0.50 

■1.50 
■1.49 
•1.43 
■1.47 
■1.46 
■1.45 
•1.44 
•1.43 
■  1    42 

■1.41 
■1.40 


-1.68 
-1.68 
-1.68 
-1.67 
-1.66 
-1.65 
-1.64 
-1.63 
-1  .62 
-1.61 
-1.60 


D 


■1.82 

■1.82 
•1.82 
•1.82 

■1.81 
■1.80 
■1.79 
•1.78 
■1 .77 
■1.76 
■1.75 


Water  Salinities    (g/kg) 


Depth 
(m) 

Profile 

A 

B 

C 

D 

0 

28.00 

28.00 

31.00 

33.51 

5 

28.50 

28.01 

31.01 

33.52 

10 

29.00 

28.02 

31.02 

33.53 

15 

29.50 

28.03 

31.06 

33.55 

20 

30.00 

28.04 

31.15 

33.57 

25 

30.50 

28.05 

31.25 

33.60 

30 

31.00 

28.06 

31.35 

33.63 

35 

31.50 

28.07 

31.45 

33.66 

40 

32.00 

23.08 

31.55 

33.69 

45 

32.50 

28.09 

31.65 

33.72 

50 

33.00 

28.10 

31.75 

33.75 
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Table  II  presents  a  summary  of  runs  conducted  with  their  respective 
variable  combinations. 

D.  MODEL  LIMITATIONS 

In  addition  to  the  assumptions  made  earlier  in  arriving  at  the 
simplified  model  of  an  open  lead,  basic  limitations  in  the  model's 
capabilities  are:  1)  the  atmosphere  is  treated  as  an  external  parameter 
such  that  it  affects  lead  dynamics  but  cannot  be  affected  in  turn  by  the 
lead;  2)  sea  ice  is  treated  as  a  quasiexternal  parameter  since  it  is 
depicted  as  a  simple,  uniform  mass  floating  on  the  water's  surface. 
Density  and  salinity  are  constant  throughout  this  mass  in  both  space  and  time, 
and  the  vertical  temperature  gradient  is  taken  to  be  linear  throughout.  Its 
thermal  conductivity  is  determined  only  at  the  ice-water  boundary,  which  is 
warmest  and  therefore  limiting  through  the  relation  given  in  equation  (43), 
and  an  upper  limit  is  applied  to  the  upward  heat- flux  density  through  the 
ice  such  that  at  no  time  is  the  flux  through  the  ice  allowed  to  exceed  that 
over  open  water;  3)  mechanical  lead  closure  (such  as  by  wind-drifted  ice) 
is  disallowed;  4)  snow  cover  over  the  lead's  surface  is  not  modelled;  5) 
super-cooling  of  surface  waters  [Coachman  1966]  is  not  portrayed;  6)  as 
convection-induced  mixing  occurs  instantaneously  in  the  model,  so  also 
does  surface  ice  melting  whenever  heat  is  brought  to  the  surface  by 
convection.  This  results  in  an  artificiality  in  that,  although  heat  and 
salt  are  conserved  in  the  model,  heat  brought  to  the  surface  under  an  ice 
cover  immediately  melts  ice  until  no  further  surface  heat  is  available  for 
melting.  Therefore,  under  this  limitation,  it  is  not  possible  for  surface 
heat  to  be  advected  downstream  under  an  ice  cover  to  promote  downstream 
melting.  Only  in  cases  when  no  ice  is  present  to  be  melted  may  this  heat 


37 


Table  II.  Schedule  of  Runs  and  Variable  Combinations 


Run 

Variable  Profile 

Period  of  Integration 
(days) 

U 

T 

S 

1 

A 

A 

A 

2 

2 

A 

A 

B 

2 

3 

A 

B 

A 

2 

4 

A 

B 

B 

2 

5 

B 

A 

A 

2 

6 

B 

A 

B 

2 

7 

B 

B 

A 

2 

8 

B 

B 

B 

2 

9 

C 

A 

B 

2 

10 

D 

A 

B 

2 

11 

E 

C 

C 

2 

12 

E 

D 

D 

2 

13 

£ 

n 

n 

30 

14 

F 

A 

B 

2 
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be  advected  downstream;  7)  in  situations  wherein  vertical  salinity  gradients 
existing  in  the  water  column  are  weaker  than  a  certain  minimum  gradient, 
the  model  will  over-mix  the  water  column  in  depth.  Specifically,  the  model 
forms  ice  in  response  to  surface  cooling,  thereby  rejecting  salt  into  the 
surface  water  layer  of  depth  Z.  The  model  then  mixes  the  induced  un- 
stable region  as  determined  by  relative  sigma-t  values  of  density.  If 
however,  there  is  heat  available  at  some  shallower  depth  than  that  to 
which  the  water  column  was  mixed  which  could  have  decreased  the  initial 
amount  of  ice  formed,  then  although  salt,  ice,  and  heat  are  ultimately 
conserved,  potential  energy  is  not  conserved.   Initial  freezing  and  resul- 
tant salt  rejection  have  induced  convection  to  a  greater  depth  than  would 
have  actually  been  the  case  in  a  continuous  process.  Therefore,  for 
vertical  salinity  gradients  weaker  than  a  certain  minimum  gradient,  the 
model's  numerical  convective  process  acts  as  a  potential  energy  source. 
Quantitatively,  the  limiting  minimum  vertical  salinity  gradient  which  can 
be  accurately  treated  by  this  model,  in  g/kg  per  vertical  grid  space,  is 

Pi  .  (Sw  -  S^  •  AI/A-Z  (45) 

where  p.  is  the  ice  density,  Sw  is  the  salinity  of  sea  water,  S..  is  the 
salinity  of  sea  ice,  A  I  is  the  incremental  change  in  the  volume  of  ice, 
and  AZ  is  the  vertical  grid  spacing.  In  this  model,  the  limiting  gradient 
is  in  the  neighborhood  of  1  X  10"3  g/kg  per  vertical  grid  space,  which  is 
lower  by  an  order  of  magnitude  than  the  weakest  salinity  gradient  used. 
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IV.  PRESENTATION  OF  DATA 

Some  results  of  integrating  the  model  using  the  combinations  of 
environmental  variables  summarized  in  sections  III . B.  and  III .C.  are 
presented  in  tabular  and  graphic  form.  Three  primary  output  variables  are 
treated:  surface  ice  thickness,  surface  heat  loss  to  the  atmosphere,  and 
depth  of  penetration  of  surface-freezing  induced  convective  overturn. 
Throughout  the  section,  these  three  output  variables  are  related  by  run 
number  (section  III.C.)  to  the  input  variables  of  current  magnitude, 
vertical  temperature  profile  and  vertical  salinity  profile  (section  III.B.). 

A.  TEST  CASES 

Ice  thicknesses  are  treated  in  two  ways:  by  time  growth  rate  and  by 
spatial  profile  across  the  lead.  Values  for  growth  rates  were  arrived  at 
by  tabulating  the  maximum  ice  thickness  occurring  anywhere  in  the  lead  as 
a  function  of. time.  Values  for  spatial  profiles  were  taken  as  the  ice 
thicknesses  across  the  model's  surface  (200  meters)  after  48  hours  of 
freezing.  Values  over  the  right  50  meters  of  the  model's  surface,  initially 
covered  by  200  cm  of  ice,  are  tabulated  to  show  the  downstream  surface 
effects  of  lead  thermodynamics  on  the  equilibrium  ice  cover. 

Surface  heat  losses  to  the  atmosphere  are  presented  as  cumulative  heat- 
loss  densities  (heat  loss  per  unit  surface  area)  occurring  at  each  grid 
point  over  the  surface  of  the  initially  open  lead  after  48  hours  of  lead 
reclosure.  Additionally,  these  surface  heat-loss  densities  are  averaged  and 
tabulated  separately  as  mean  surface  heat-loss  densities. 
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Finally,  maximum  depths  of  penetration  of  surface-freezing  induced 
convective  overturn  are  tabulated.  Due  to  the  numerical  nature  of  the 
model  however,  these  depths  should  be  treated  as  bounding  depths  rather 
than  precise  depths  as  would  be  indicated  for  continuous  schemes. 

B.  REAL  CASES 

Ice  thicknesses,  surface  heat  losses  to  the  atmosphere,  and  depths  of 
penetration  of  surface-freezing  induced  convective  overturn  are  treated  in 
an  analogous  manner.  Additionally  however,  ice  thickness  and  mean  surface 
heat-loss  densities  are  presented  as  functions  of  time  for  run  number  13, 
in  which  the  model  was  integrated  to  represent  30  days  of  lead  reclosure. 
Values  for  surface  ice  thicknesses  were  determined  analogously  to  those  of 
48-hour  ice-growth  rates. 
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Table  VI.  Mean  Surface  Heat-Loss  Densities,  Test  Cases, 


48-Hour  Mean-Surface  Heat  Losses 
(x  102  cal/cm2) 


Run 

Number 

A0 

1 

14.9 

2 

17.0 

3 

13.7 

4 

14.9 

5 

14.9 

6 

23.2 

7 

13.7 

8 

14.5 

Table  VII.  Maximum  Depth  of  Penetration  of  Convective  Overturn, 
Test  Cases. 


Run 

Depth 

Number 

(m) 

1 

5   • 

2 

10 

3 

5 

4 

10 

5 

5 

6 

20 

7 

5 

8 

20 

51 


CO 
cd 

CO 

03 
O 


to 

CD 


00 

cd 
+-> 
03 
q: 

x: 
O 

s_ 

CO 

cd 
u 


CO 

r— 
fO 

> 

CD 


i- 

ZJ 

o 


en 


+-> 

«3 

<u 
o 

03 

<+- 

S- 
15 
CO 

CO 

-a 

rt> 

a; 

_i 

CD 
4-> 


CO 

c 

•• — 

i- 

s- 

O 

u 
o 


CO 
co 

CD 

C 

u 


CD 

o 


X 

03 


CD 
03 


r-^ 

MD 

00 

• 

• 

•vT 

«3- 

«3- 

00 

co 

CO 

C\J 

O 

• 

• 

"=d" 

OJ 

CM 

00 

CM 

r^ 

CO 

c 

CNJ 

• 

• 

o 

CO 

CD 

CT> 

•  ^ 

r— 

h~ 

4-> 

03 

CD 

c 

CD 

CD 

~o 

|-» 

CO 

03 

«3- 

• 

• 

CD 

C\l 

CO 

CO 

_J 

•"" 

■— 

S- 

CD 
+-> 

<+- 

< 

CD 

*— 

evi 

r~ 

•  i — 

<x> 

• 

. 

t— 

1 — 

CO 

CO 

st 

00 

co 

C\J 

00 

CO 

0 

0 

0 

S- 

CD 

c  -a 

1 — 

C\J 

Z3    E 

1 — 

r— 

o;  3 

2: 

52 


CO 


o 
*3- 


C\J 
CO 


<3- 


c 
o 


03 

s_ 
<v 

c 

CD 

-a 

<D 


CD 

+-> 

4- 

< 


UD 


-        CO 


C\J 


o 

CM 


in 


ID 


(LUD)     SSSU>jDLlJl    301    WniUlXPl 


c;-5 


Table  IX.  Ice  Thickness  Profiles,  Real  Cases 


Ice  Thicknesses  (cm)  Across  Lead  Surface 
48  Hours  After  Lead  Generation. 


X 

(m) 

Run  Number 

11 

12 

10 

24.7 

24.6 

20 

24.7 

24.6 

30 

24.7 

24.6 

40 

24.7 

24.6 

50 

24.6 

24.3 

60 

24.3 

24.0 

70 

24.4 

24.2 

80 

24.4 

24.2 

90 

24.4 

24.1 

100 

24.4 

23.9 

110 

24.3 

23.9 

120 

24.2 

23.6 

130 

24.0 

23.5 

140 

24.0 

23.4 

150 

203. 

203. 

160 

203. 

203. 

170 

203/ 

203. 

180 

203. 

203. 

190 

203. 

203. 

200 

203. 

203. 
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Table  X.     Surface  Heat-Loss  Densities,   Real   Cases 

48-Hour  Surface  Heat  Losses 
(x  102  cal/cm2) 


X 
(m) 

Run  Number 

11 

12 

10 

13.6 

13.9 

20 

13.6 

13.9 

30 

13.6 

13.9 

40 

13.6 

13.9 

50 

13.7 

14.3 

60 

14.2 

14.7 

70 

14.0 

14.4 

80 

13.9 

14.3 

90 

13.8 

14.3 

100 

13.8 

14.3 

110 

13.8 

14.3 

120 

13.8 

14.3 

130 

13.8 

14.3 

140 

13.8 

14.2 

150 

1.99 

2.01 

160 

1.99 

2.01 

170 

1.99 

2.01 

180 

1.99 

2.01 

190 

1.99 

2.01 

200 

1.99 

2.01 
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Table  XI.  Mean  Surface  Heat-Loss  Densities,  Real  Cases. 


48-Hour  Mean  Surface  Heat  Losses  (x  102  cal/cm2) 


Run 
Number 

aQ 

11 
12 

13.3 
14.2 

Table  XII.  Maximum  Depth  of  Penetration  of  Convective  Overturn, 
Real  Cases. 


Run 

Number 

Depth 
(m.) 

11 
12 

15 
15 
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Table  XIII.     Ice  Growth  Rate,  Run  13. 


Maximum  Ice  Thickness   (cm)  Occurring  in  the 
Lead's  Surface  at  Five-Day  Intervals. 


Ice 

Day 

'  Thickness 

5 

40.6 

10 

59.6 

15 

74.9 

20 

88.3 

25 

101. 

30 

112. 

Table  XIV.     Mean  Surface  Heat-Loss  Densities,   Run  13. 

9  9 

Mean  Surface  Heat-Loss   Densities   (xlO^  cal/cnf) 
at  Five-Day  Intervals. 


Day 

aQ 

5 

23.5 

10 

34.7 

15 

43.7 

20 

51.6 

25 

58.8 

30 

65.5 
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V.  CONCLUSIONS 

Through  the  preceding  data  and  the  time-dependent  temperature-salinity- 
density  fields  generated  in  the  lead  model,  numerous  insights  to  the  process 
of  lead  reclosure  were  afforded,  as  well  as  to  the  hydrodynamic  validity  of 
the  model  under  certain  conditions. 

Perhaps  most  pronounced  was  the  very   strong  role  played  by  the  near- 
surface  vertical  salinity  gradient  in  maintaining  water-column  stability. 
A  strong  halocline  near  the  surface  tended  to  insulate  deeper  waters  in  the 
lead  from  the  thermodynamic  activity  occurring  at  the  surface.  Surface 
waters  whose  densities  had  been  increased  by  salt  rejected  by  newly  forming 
ice,  was  not  able  to  mix  vertically  until  density  instabilities  occurred  in 
the  surface  region,  thus  inducing  convective  overturn.  Since  the  presence 
of  a  strong  surface  halocline  would  require  higher  surface-salinity  accumula- 
tions and  hence  permit  more  ice  formation  and  associated  salt  rejection 
before  this  surface  layer  would  become  unstable,  a  strong  halocline  would 
be  expected  to  inhibit  convective  overturn,  both  in  frequency  of  occurrence 
and  depth  of  penetration  in  the  lead.  Added  to  this  insulating  effect  is 
the  fact  that,  as  surface  salinities  increase,  surface  freezing  temperatures 
correspondingly  decrease  such  that  freezing  rates  are  proportionately  lowered. 

The  role  played  by  the  vertical  temperature  gradient  in  lead  thermodynamics 
was  strongly  dependent  upon  the  existing  vertical  salinity  gradient.  Since 
density  changes  in  sea  water  at  temperatures  near  the  freezing  point  are 
almost  totally  a  function  of  salinity  changes,  the  insulating  qualities  of 
the  thermocline  are  undefinable  when  treated  independently  of  the  halocline. 


However,  in  situations  wherein  weak  vertical  salinity  gradients  (poor  in- 
sulation) and  strong  vertical  temperature  gradients  coexisted,  convective 
overturn  resulted  in  relatively  large  quantities  of  heat  being  transported 
to  the  surface.  Freezing  and  associated  salt  rejection  were  thus  inhibited 
and  surface  melting  was  induced.   In  this  manner,  the  temperature  stratifi- 
cation was  superimposed  on  the  primary  control  exerted  by  the  vertical 
salinity  gradient. 

The  effect  of  horizontal  current  magnitude  on  lead  stability  appeared, 
as  in  the  case  of  the  thermocline,  to  be  secondary  to  the  effect  of  the 
halocline.  Strong  currents  will  advect  surface  salt  accumulations  out  of 
the  lead  relatively  quickly,  and  will  in  general  allow  less  time  for  lead 
thermodynamics  to  affect  a  given  parcel  of  water  advectively  transiting 
the  lead's  width.  This  allows  less  change  in  the  parcel's  physical  charac- 
teristics. Through  these  effects,  the  current  also  acts  as  an  insulating 
agent. 

Ice  growth  rates  were  clearly  dependent  upon  combinations  of  the  above 
three  variable  parameters.  Neglecting  spacial  ice  profiles  across  the  lead 
temporarily  and  examining  only  the  maximum  ice  thicknesses  occurring  on  the 
lead's  surface  as  a  function  of  time,  the  highest  growth  rates  were  exhibited 
for  cases  where  a  strong  current,  a  weak  thermocline,  and  a  weak  halocline 
occurred  together.  Conversely,  lowest  growth  rates  resulted  from  the  co- 
existence of  a  strong  current,  a  strong  thermocline,  and  a  strong  halocline. 
Very   little  difference  however,  was  noted  in  the  high  growth-rate  case 
when  currents  were  weak  instead  of  strong,  nor  in  the  low  growth-rate  case 
when  currents  were  weak  or  when  currents  were  weak  and  the  halocline  was 
weak.  Analysis  of  these  results  showed  that  high  growth  rates  were  in  all 
cases  associated  with  v/eak  vertical  temperature  gradient:   and  low  growth 
rates  with  strong  vertical  temperature  gradients.  Further,  that  current 


•  63 


strength  was  relatively  unimportant  in  changing  growth  rates  when  either 
strong  haloclines  were  present,  or  when  weak  temperature  gradients  together 
with  weak  salinity  gradients  were  present. 

The  effect  of  the  near-surface  vertical  temperature  gradient  appears 
then  to  dominate  in  determining  surface  freezing  rates.  Strong  gradients 
resulted  in  lower  rates  of  ice  formation  than  in  the  case  of  weak  gradients. 
Vertical  transport  of  heat  from  depth  in  the  lead  by  both  eddy  diffusion  and 
convective  overturn  is  then  modified  by  both  current  strength  and  the  near- 
surface  halocline  acting  as  insulating  agents. 

Ice  thickness  profiles  across  the  lead's  width  varied  considerably  with 
different  combinations  of  the  three  basic  variable  parameters.  Uniform  or 
nearly-uniform  ice  thicknesses  across  the  lead  were  in  all  cases  associated 
with  strong  near-surface  haloclines.  Variations  occurred  where  weak  halo- 
clines existed,  growing  larger  when  weak  haloclines  and  strong  thermoclines 
occurred  together.  The  greatest  surface  irregularity,  both  by  amplitude  of 
variation  as  well  as  by  spacial  frequency  of  variation,  was  found  when  weak 
currents,  weak  haloclines,  and  strong  thermoclines  occurred  together. 

Analysis  of  these  cases  indicated  that  the  halocline  was  the  dominant 
variable  in  determining  uniformity  of  ice  thickness  across  the  lead  through 
insulation  of  the  surface  from  the  heat  at  depth  in  the  lead.  When  convective 
overturn  occurred  as  a  result  of  increases  in  surface  salinity,  surface  freezing 
rates  were  altered  at  the  point  of  convective  overturn  and  produced  surface 
irregularities.  Since  rejected  salt  is  being  continuously  advected  down- 
stream, there  will  exist  horizontal  salinity  gradients  which,  when  exceeding 
the  magnitudes  of  vertical  salinity  gradients,  will  induce  convective  overturn 
and  resultant  surface  melting  at  that  point.  Downstream  of  this  point, 
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surface  waters  are  diluted  due  to  the  upstream  melting,  thus  experiencing 
increased  freezing  rates  since  their  freezing  point  has  effectively  been 

raised. 

The  near  surface  thermocline  appeared  to  affect  surface-ice  uniformity 
through  determining  the  amount  of  heat  transported  to  the  surface  during 
convective  overturn,  and  thus  the  amount  of  melting  and  dilution  occurring. 
Finally,  current  strength  appeared  to  act  analogously  to  the  halocline  as  an 
insulating  agent,  though  less  effectively.  By  flushing  surface-salt  accumu- 
lations out  of  the  lead  more  rapidly,  a  strong  current  appears  to  inhibit 
convective  overturn  and  hence  variations  in  ice  uniformity.  By  decreasing 
the  amount  of  change  a  given  parcel  of  water  experiences  over  a  given  hori- 
zontal distance  in  the  lead,  a  strong  current  appears  to  produce  less 
intensive  but  more  extensive  lee  effects  from  surface  dynamics. 

Surface  heat-loss  densities  varied  across  the  lead's  surface  in  response 
to  ice  thickness  profiles.  Since  open-water  sensible,  latent,  and  radiative 
heat-flux  densities  used  in  this  work  totalled  2.00  x  10"  cal/cm2-sec  as 
compared  to  a  corresponding  value  over  the  200  cm  of  ice  of  1 .15  x  10" 
cal/cm2-sec,  resultant  heat-loss  densities  were  considerably  larger  over 
regions  of  the  lead  surface  covered  with  relatively  thinner  ice.  Uniformity 
of  surface  heat-loss  densities  over  the  lead's  surface  may  then  be  directly 
related  to  the  uniformity  of  the  ice  profile  over  the  surface,  while  relative 
magnitudes  of  heat-loss  densities  may  be  inversely  related  to  relative  ice 

thicknesses. 

The  extreme  case  of  lead  instability  and  non-uniformity  (weak  current, 
strong  thermocline,  and  weak  halocline)  requires  special  comment.  A  high 

mean  surface  heat-loss  density  characterized  this  case,  which  was  expected  as 

a  result  of  its  highly  non-uniform  surface-ice  profile.  However,  the  lack 
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of  evidence  of  these  highly  non-uniform  ice  covers  occurring  in  nature  in  the 
absence  of  mechanical  deformation  forces  prompted  further  analysis.  By 
studying  model  output  fields  of  runs  2,  6,  9,  10,  and  14  (particularly  the 
horizontal  density  gradients  generated  in  the  lead),  an  upstream  pressure 
force  was  found  to  exist  in  the  vicinity  of  the  source  boundary  which  was 
capable  of  performing  on  the  order  of  ten  ergs  of  work  against  a  unit  mass 
of  water  entering  the  lead  at  the  source  boundary.  Since  a  thermodynamically 
generated  force  of  this  type  would  completely  alter  the  hydrodynamics 
occurring  in  the  lead  for  low  current  magnitudes,  the  model  appears  to  be  non- 
representative  of  actual  lead  thermodynamics  when  using  very  low  current 
magnitudes  with  strong  thermoclines  and  weak  haloclines.  Stronger  currents, 
on  the  other  hand,  had  generally  weaker  horizontal  density  gradients  associated 
with  them  and  were  also  proportionately  less  affected  by  these  pressure  forces 
since  their  capacity  to  do  work  against  the  horizontal  pressure  force  in- 
creases as  the  square  of  their  magnitude. 

Results  of  run  14  proved  to  be  of  particular  interest  through  the  in- 
sights they  afforded  to  the  process  of  convective  overturn.  Since  no  current 
was  present,  there  was  no  horizontal  coupling  between  model  grid  points.  As 
a  result,  although  freezing  progressed  uniformly  over  the  surface  of  the  lead, 
strong  differential  horizontal  density  gradients  developed  along  the  'source' 
boundary  as  cooling  progressed.  Of  prime  interest  however,  was  the  graphic 
illustration  of  the  role  played  by  the  current  as  a  convective-overturn 
scaling  factor.  With  zero  current,  a  strong  temperature  gradient  and  a 
weak  salinity  gradient,  freezing-induced  convective  overturn  penetrated  to 
a  depth  of  35  meters  roughly  18  hours  after  lead  generation.  Without  the 
flushing  action  of  a  horizontal  current  or  the  insulating  effect  of  a  strong 
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halocline,  the  effects  of  surface  thermodynamics  were  felt  much  deeper  in 
the  lead.  It  appears  probable  then,  that  in  cases  wherein  weak  currents  and 
weak  haloclines  exist  in  open  leads  (marginal  vertical  stability),  surface 
freezing  will  induce,  through  convective  overturn  and  vertical  mixing,  sub- 
stantial vertical  and  horizontal  gradients  in  both  temperature  and  salinity. 
As  a  result,  significant  fluxes  of  heat,  salt,  and  momentum  would  be  genera- 
ted in  the  waters  of  the  lead,  and  sophistication  beyond  that  possessed  by 
the  work  at  hand  would  be  required  to  accurately  represent  the  resultant 

lead  hydrodynamics. 

Summarizing  the  analysis  of  results  from  integrating  the  open-lead  model 
over  the  experimental  ranges  of  variables,  the  near-surface  vertical  salinity 
gradient  appears  to  exert  the  dominant  force  in  maintaining  vertical  density 
stability  in  the  water  column  through  its  insulating  effect.  By  inhibiting 
convective  overturn,  a  strong  halocline  was  associated  with  horizontally- 
uniform  ice  thickness  profiles  and  consequently  uniform  surface  heat-flux 
densities.  Current  strength  and  the  vertical  temperature  gradient  appeared 
to  be  secondary  as  insulating  agents. 

The  near-surface  temperature  gradient  appeared  to  be  the  dominant 
variable  in  controlling  the  rate  of  ice  formation  since  its  strength 
directly  influences  the  amount  of  heat  transported  to  the  ice-water  boun- 
dary from  depth  by  both  eddy  diffusion  and  convective  overturn.  In  this 
case,  the  halocline  and  the  current  strength  appeared  secondary  in  their 
effects  on  the  freezing  rate. 

Mean  surface  heat-loss  densities  were  higher  for  more  non-uniform 
surface  ice  thicknesses,  which  in  turn  may  be  inversely  related  to  the 
strength  of  the  near-surface  halocline.  This  conclusion  could  prove  to  be 
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of  significance  to  critical  examinations  of  the  heat  budget  over  large 
areas  of  the  arctic  in  which  significant  regional  variations  in  the  near- 
surface  halocline  might  reasonably  be  expected. 

Representative  freezing  rates,  ice-thickness  profiles,  and  surface 
heat-loss  densities  found  in  integrating  the  model  using  real-case  variables 
agreed  well  with  values  in  the  literature.  Due  to  the  similarity  of  current 
strengths,  near-surface  thermoclines,  and  near-surface  haloclines,  the  fore- 
going relations  were  generally  not  discernable,  although  the  effect  of  down- 
stream salt  accumulation  was  noticed  in  slightly  decreasing  downstream  ice 
thicknesses. 

Finally,  it  is  concluded  that,  with  the  exception  of  cases  in  which 
weak  horizontal  currents  occur  together  with  strong  thermoclines  and  weak 
haloclines,  the  model  presented  herein  provides  an  accurate  physical 
representation  of  the  thermodynamics  occurring  in  an  open  lead  in  central 
arctic  sea  ice  over  the  experimental  range  of  variables. 
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VI.     RECOMMENDATIONS   FOR  FURTHER  RESEARCH 

Although  the  surface  boundary  of  the  model   could  be  considerably  more 
sophisticated  in  its  treatment  of  energy  fluxes  and  transport  processes, 
models  already  exist  which  have  this  sophistication,  notably  that  of  Maykut 
and  Untersteiner  [1969,  1971].     Coupling  the  model   at  hand  to  their  model 
would  therefore  seem  a  more  reasonable  approach  to  refinement  of  the  surface 
boundary  in  order  to  model   air-ice-sea  interactions  and  mass-energy  balances. 

The  hydrodynamics  occurring  in  a  refreezing  lead  however,  are  still   not 
accurately  modelled,  nor  are  they  fully  understood  despite  the  insights 
afforded  by  the  research  performed.      Extended  development  of  this  work 
would  promote  further  conceptual   understanding  of  these  processes  as  well 
as  provide  a  more  physically-accurate  model.      Incorporation  of  accelerations 
resulting  from  thermodynamically-generated  horizontal   pressure  gradients, 
as  well   as  the  coriolis  acceleration,  would  considerably  refine  the  model. 
Additional    refinement  would  be  gained  by  removing  previously  discussed 
limitations  on:     1)  the  minimum  vertical    salinity  gradient  able  to  be 
accurately  treated,  and  2)  the  preclusion  of  downstream  advection  of  warmer- 
then-freezing-point  waters  under  an  ice  cover. 

In  summary  then,  further  research  could  profitably  be  applied  to  the 
work  at  hand  to  further  refine  the  hydrodynamic  accuracy  of  the  model,  as 
well   as  to  couple  it  to  larger-scale  models  of  air-ice-sea  interactions  and 
mass-energy  balance  predictive  schemes. 
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APPENDIX  A 
COMPUTER  PROGRAM  DESCRIPTION 

The  computer  program  is  written  in  FORTRAN  IV  and  was  used  with  the 
IBM  360/67  computer  system  at  the  W.R.  Church  Computer  Center,  Naval 
Postgraduate  School.  The  program  has  been  generalized  such  that  sea 
water  temperature  and  salinity  profiles  used  as  initial  and  left  boundary 
fields  are  treated  as  input  data.  A  design  feature  of  the  program  allows 
the  use  of  an  internally  variable  time  step  to  facilitate  adjusting  the 
model's  temporal  resolution.  Additionally,  the  model  allows  variable  grid 
spacing  as  well  as  variable  lead  dimensions,  requiring  changes  only  in  the 
DIMENSION,  COMMON  and  FORMAT  statements.  Although  all  currents  used  in 
the  model  were  constant  in  time  and  space,  the  model  allows  use  of  both  a 
depth-dependent  and  time-dependent  horizontal  current,  again  internally 
variable.  With  only  minor  modifications,  the  program  could  be  easily  ex- 
tended to  include  a  bottom-boundary  heat  and  salt  flux,  more  sophisticated 
treatment  of  air-sea  interactions  at  the  surface  boundary,  and  feedback  of 
accelerations,  resulting  from  horizontal  density  gradients  generated  by  lead 
thermodynamics,  to  the  motion  of  water  through  the  lead. 

The  program  is  divided  into  two  basic  sections:  the  main  program  and 
subroutine  OVRTRN.  The  main  program  is  further  subdivided  into  seven  parts: 
problem  initialization,  data  input,  computation  of  time-dependent  temperature 
and  salinity  fields,  ice-accretion/melting  and  salt  rejection/dilution,  with 
vertical  stability  check  and  computation  of  resultant  heat  flux  densities, 
filling  temperature  and  salinity  arrays  with  resultant  values,  computation 
of  applied  surface  heat  flux  densities  and  conversion  to  an  integrated  surface 
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temperature  gradient,  and  writing  output.  Subroutine  OVRTRN  accomplishes 
vertical  convective  overturn  and  mixing  of  unstable  regions  in  the  water 

column. 

FORTRAN  IV  symbols  for  primary  program  parameters  are  defined  below: 
AH  Downstream-advected  heat-flux  density  (cal/cm  -sec) 

AHF        Applied  surface  heat-flux  density  (cal/cm  -sec) 

BTC        Sea  ice  thermal   conductivity  constant  (cal-cm  /g-sec) 

2 
CAH         Cumulative  downstream-advected  heat  density  (cal/cm  ) 

CAHF  Cumulative  applied  surface  heat  density   (cal/cm  ) 

CHC  Cumulative  subsurface  cooling  heat  density   (cal/cm  ) 

CHF  Cumulative  applied  surface  heat  density  less  cumulative 

downstream-advected  heat  density   (cal/cm2) 

CHIF  Cumulative  heat  density  of  surface-ice  formation   (cal/cm  ) 

CI  Thermal   conductivity  of  sea  ice   (cal/cm-sec-°C) 

CL  Sea-water  chlorinity   (g/kg) 

CO  Thermal    conductivity  of  pure  ice   (cal/cm-sec-°C) 

CP  Specific  heat  of  sea  water  (cal/g-°C) 

DI  Rate  of  surface  ice  formation    (cm3/cm2-sec) 

DSALT  Salt  rejected  per  cm3  of  ice  formed   (g/kg) 

DTAU  Time  step   (sec) 

DTEMP  Air-water  temperature  difference   (°C) 
DX  Horizontal   grid  spacing   (cm) 

DZ  Vertical   grid  spacing   (cm) 

EPS        Allowable  rate-of-freezing  error   (cnrVcm  -sec) 
ERR        Allowable  ice  thickness  error  (cm  /enr) 
GRADI     Vertical    ice  temperature  gradient  (°C/cm) 
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H  Model   depth  extent   (cm) 

2 


HC 


Subsurface  cooling  heat-flux  density  (cal/cm  -sec) 
HF  Applied  surface  heat-flux  density  less  downstream-advected 

heat  flux  density   (cal/cm  -sec) 
HIF        Heat-flux  density  of  surface-ice  formation   (cal/cm  -sec) 
HO  Heat  of  formation  of  pure  ice   (cal/g) 

HS  Heat  of  formation  of  sea  ice   (cal/g) 

MM         Number  of  grid  columns 
MV1         Number  of  grid  rows 
MVH1       Total   number  of  grid  points 
QI  Heat-flux  density  through  sea  ice  to  the  atmosphere 

(cal/cnr-sec) 
QL  Open-water  latent  heat-flux  density  to  the  atmosphere 

(cal/cnr-sec) 

QR  Open-water  radiative  heat-flux  density  to  the  atmosphere 

(cal/cnr-sec) 
QS  Open-water  sensible  heat-flux  density  to  the  atmosphere 

(cal/cm  -sec) 

RHO  Density  of  sea  water  (g/cm  ) 

RHOI  Density  of  sea  ice   (g/cm3) 

S  Sea-water  salinity   (g/kg) 

SALI  Sea-ice  salinity   (g/kg) 

SI  Sea-ice  salinity  (g/cm3) 

SIGT  Sigma-t 

SSCHG  Salinity  change  in  surface  layer  of  depth  DZ   (g/kg) 
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SSFLUX     Salinity  change  in  one  cm3  of  sea  water  at  the  ice-water 

boundary   (g/kg) 
STGD        Total   applied  surface  temperature  gradient  at  the  water's 

surface  (°C/cm) 
T  Sea-water  temperature   (°C) 

Tau  Time   (sec) 

TF  Sea-water  freezing  temperature   (°C) 

TGRADL     Surface  temperature  gradient  due  to  latent  heat-flux 

density  (°C/cm) 
TGRADR     Surface  temperature  gradient  due  to  radiative  heat- flux 

density   (°C/cm) 
TGRADS     Surface  temperature  gradient  due  to  sensible  heat-flux 

density   (oc/cm) 
TI  Temperature  of  lower  ice  surface   (°C) 

U  Horizontal   current  magnitude   (cm/sec) 

VEC  Vertical   eddy  conductivity   (cal/cm-sec-°C) 

VSK  Vertical    salt  eddy  diffusivity   (cnf/sec) 

VTK  Vertical    temperature  eddy  diffusivity   (cm  /sec) 

W  Lead  width   (cm) 

Z  Depth    (cm) 

11  Surface  ice  thickness   (cm3/cm2) 

A  descriptive  flow  diagram  of  the  program  is  shown  in  Figure  16. 
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Figure  16.  Descriptive  Flow  Diagram  of  Program. 


Problem  Initialization: 
DIMENSION/COMMON 
Function  statements, 
Constants 


Input  of  Initial  Conditions: 
Temperature  TNEW, 
Salinity  SNEW, 
Sigma-t  SIGT, 
Equilibrium  Ice  Thickness  ZI 


Computation  of  Time-Dependent 
Temperature  and  Salinity  Fields 
TNEW,  SNEW 


Computation  of  Surface  Temperature  Chances 
THF 


(202)—  —  —  — 
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Figure  16   (continued) 


202 
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Computation  of: 
Surface  freezing  temperatures  TF, 
Surface  temperature  changes  due  to  ice 
formation  TIF, 

Iterative  rate  of  ice  formation  HICE, 
Iterative  volume  of  ice  formed  ZICE, 
Iterative  salt  rejected  SSCHG, 
Resultant  surface  temperature  and  salinity 
fields  (if  melting  occurred  as  a  result  of 
convective  overturn,  the  resultant  dilution 
is  distributed  over  the  mixed  region  of  the 
water  column)  TNEW,  SNEW 


Computation  of  Surface  Ice  Thicknesses  ZI 
(iteration  to  account  for  freezing-point 
changes  as  a  result  of  salinity  changes) 


1 


Vertical  Stability  Check 


OVRTRN 


Check  for  Surface  Heat  Brought 
From  Depth  Which  Would  Promote 
Melting 


J 


Computation  of: 

Rate  of  ice  formation, 
Water-column  heat-flux  densities, 
Cumulative  water  column  heat  densities 


7k 


1 


Fill   Temperature  and  Salinity  Arrays 
T 
S 


Confutation  of  Surface  Heat-Flux  Densities  and 
Conversion  to  an   Integrated  Surface  Temperature 
Gradient 


'-i+inn     H, 


©^ 
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CENTRAL  ARCTIC  OPEN-LEAD  MODEL*  ,„  , 

DIMENSION  T(232)  ,S(232)  ,PUF<22)  fAH(22) ,CAH(22) , 
1H  F(22,CH  H22   DK22I  Z I < 22 ) , TF < 22 ) , HF < 22  )  , CHF( 22 ) , 
2STGD(22),HC(22),CHC(22>,AHF(22),CAHF(22),HUF(22) 
CCMMON  TMEW(232) ,SNEW(232) , SIGT<232) 

AT1(bI=1.-DTAU*(U(B)/DX+2.*VTK/DZ**2) 
AS1(B)=1.-DTAU*(U(B)/DX+2.*VSK/DZ**2) 

AT2(B)=OTAU*U( B)/DX 
AS2(5)=OTAU*U( e )/DX 
AT3(D)=D*VTK/DZ**2 

AS3(D)=D*VSK/DZ^"2 

SOu7  =  -0^5&9  +  lo4708*A-O.OCa57G*A**2+OcG000398*A**3 
AT (A    =A*(4o78O7-G.G98135*A+Jc0C10843*A**2)*9o001 
BT(A)=A*( 18.U3O0.  8164*A+0.01667*A**2)*J. OOOOOl 
ST(A)=-(  A-3o^8  )**2/5C3.570*(  A-»-283.)/(A  +  67e26) 
SGTU)=A+(SIG0+O.1324)*(1.0-ATMP+(i}TMP*<SIG0-0.1324) ) ) 
SH( A )= 1.005-0. 0041 3 &*A+C. 000 l098*A**2-0. 000001324* A**3 

DE(A)=1.0+< A/1000. ) 

F<  A  )=0.G-O.G966*A+0.  0000052*  A**3) 

R(C)=DSALT*C 

ERR(D)=D*EPS 

BTC=0.28 

C0-0.0C486 

DTAU=9C0 

DTEMP=25. 

DX=1000.0 

DZ=5C0.0 

EPS=1.0E-8 

H=5C"0.0 

H0=79.77 

QL=3oOE-3 

QR=2.0E-3 

QS  =  1.5E-2 

RH0I=0.91 

A  L 1  —  o  .  v 
SI=RHOI*SALI/1000.0 
VSK=10oO 
VTK=10.0 
W= 20000.0 
MH=IFIX(W/OX) 
MV=IFI X(H/DZ) 
MV1=MV+1 
KF1=MH+1 
MVH1=MV1*MH1 

1  READ( 5t 100JIRUN 

IF( IRUN.EQ.99)G0    TO    1000 

READ (5,  15';)  (TNEWd  )  ,I  =  1,MV1) 

R E  AD  (  5  ,  1  50  )  (  SN  E  vi  (  I  )  » I  =  1 » M V 1 ) 

CL=C(SNEW(1 ) ) 

TF(  1)=F(CL) 

TNEW(1)=TF(1) 

M=MV1 

DO    3    J=2,MH1 

DO    2    K=ltMVl 

I=M+K 

TNEW( I )=TNEW(K) 

T(  I  )=TMEW(K) 

SNEWd  )=SNEW(K) 

S(  I  )=SNEW(K) 

2  CONTINUE 
M=M+MV1 

3  CONTINUE 

DO  4  1=1, MV1 
CL=C(SNEW( I) ) 
SIGO  =  Sr  (CD 
ATMP=AT(TNEW( I ) ) 
BTMP=BT (TNEW( I ) ) 
SMT=ST(TNEW( I) ) 
SIGT( I l=SGT(SMTJ 
T( I)=TNEW( I) 
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S( I)=SNEW( I) 

4  CONTINUE 

DO  5  I=12MH1 
AH< I )=O.C 
AHF(I)=0«0 
CAH(  I)=0.0 
CAHF( I )=C„0 
CHC(I)=0.0 
CHF( I)=0.0 
CHIF( I  )=2.0 
DI (I )=uoO 
HC( I)=0.0 
HF(  I)=OoO 
HIFU  )=0o0 
ZI( I )=0.0 

5  CONTINUE 
J=MHl-(MH/4) 
DO  6  I=J,MH1 
ZI(I)=20CoO 

6  CONTINUE 
TAU=O*0 
GO  TO  201 

101  TAU=TAU+DTAU 
DO  10  2  I=2»MH1 
AHF( I)=HUF( I ) 

CAHF(  I)=CAHF(  I  )  +  PUF(  I) 

102  CONTINUE 
J=MV1+1 

DO  103  1=2, MH 

Z=O.C 

Jl= J— MV1 

TNEW(  J)=AT1(Z)*TU ) +AT2C Z ) *T( Jl ) 
1+AT3(DTAU)*2.*(T( J+1)-DZ*STGD( I  )  J 

SNEW(J)=AS1(Z)*S(J)+AS2(Z)*S(J1) 
l+AS3(DTAU)*2o*S(J+l) 

Z  =  H 

O  — U  T|-|  V 

J1=J-MV1 
TNEW(J)=AT1(Z)*T(J)+AT2(Z)*TU1) 

1+AT3(DTAU)*2.*T( J- 1 ) 

SNEW( J)=AS1(Z)*S(J  )+AS2(Z)*S(Jl) 
l+AS3(DTAU)*2o*S( J-l > 

J  =  J  +  1 

103  CONTINUE 
Z=0.0 
K=MVH1-MV 

TNEw'(K)  =  ATl<Z)*T(K)-i-AT2(Z)*T(J}  +  AT3(DTAU}*2.* 

1SNEW(K)=ASl(Z?*s'(K)+AS2(Z)^S(J)+AS3(DTAU)*2,*S(K+l) 

Z=H 
K=MVH1 

J—  |/_  V)W  1 

TNEW'(K)=AT1(Z)*T(K)+AT2(Z)*T(J)+AT3(DTAU)*2«*T(K-1  j 
SNEW(K)=ASl(Z)*S(K)+AS2(Z)*S(J)+AS3(DTAU)*2e*S(K-l) 

DO  105  J=2,MH 
Z  =  DZ 

Kl=( ( J-l  )*MVl)+2 
K2=( J*MV1)-1 

DO  104  K=K1,K2 

TNEw7k)=AT1(Z)*T(K)+AT2< Z)*T (Jl )+AT3(DTAU)* 

1SNFwTk)=ASi!z)-S(K)+AS2(Z)*S(J1)+AS3(DTAU)* 
1(S(K-1)+S(K+1) ) 
Z=Z+DZ 

104  CONTINUE 

105  CONTINUE 
K1=1+(MVH1-MV) 
K2=MVH1-1 

Z  =  DZ 

DO  106  K=K1,K2 
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JrrK-MV1 

TNEW(Ki  =  ATKZ)*T(K)  +  AT2(Z)*T(J)  +  AT3(DTAU)* 

1SNEw7k)=ASi!z)*S(K)+AS2(Z)*S(J)+AS3(DTAU)^ 
1(S(K-1 )+S(K+l> ) 

Z=Z+DZ 
1C6  CONTINUE 
2C1  J=l 

DO  213  I=2,MH1 

J=J+MV1 

J3=J+MV 

MELT=G 

HIF(I)=OoC  ,  i%1 

THF=DZ/2o*(T(J  )-TNEW(J)  ) 
2C2  CL=C(SNEW( J) ) 

TF( I)=F(CL) 

TTEKP=TNEW(J) 

IF(TTEMPoLEoTF(I ))TTEMP=TF(I ) 

I^TTD^l:Si:^i!Tl^:fM«I..Lt.ERR!DTAU.».T!F=C,0 

CP=SH(SNEM( J) ) 
SIGOSC-(CL) 
ATMP=AT(TTEMP) 
BTMP=BT(TTEMP) 

SMT=ST(TTEMP) 
SIGT( J)=SGT(SMT) 
RHO=DE( SIGT( J) ) 
ZKW=RHO*CP/DTAU 
SRATIO=SALI/SNEW(J) 
HS=HC*( l.O-SRATIO) 
HICE=(TIF*ZKW)/(RHOI*HS) 
GO  TO  204 

203  HICE=-(ZI( I )/DTAU) 
TIF*=HICE*RhOI*HS/ZKW 

204  Gl=-(DTAUvHlCEJ 
IF(Gl.GT.ZKI)  >G0  TO  203 
DSALT-RHOI*(SNEW(J 5-SALI ) 
ZICE=DTAU*HICE 
SSFLUX=R(ZICE) 
SSCHG=SSFLUX/DZ 
IF(MELT.E0«C)G0  TO  2  06 
SDIST=SSCHG/JSA 

Kl=J+( JSA-1) 

DO  205  K=JtKl 

SNEW(K)=SNEW(K)+SDIST 

205  CONTINUE 
GO  TO  207 

206  SNEW( J)=SNEW( J J+SSCHG 

2C7  TNE  W ( J )  =  TNE  W( J )  +  ( T I F*2 . /DZ ) 

ZI( I)=ZI (I )+DTAU*HICE  ,  ,.,,,,,.  ,  r 

IF(ZI(  I  )  el_ToERR(  DTAO.AND.  ( HICE.LE.O.o  )  )  Z H  I  )=>./* 0 

HIF( IJ=HIF( I )+TIF*ZKW 

IF(ABS(HICE)oGTcEPS)GO  TO  2-J2 

CL=C(SNEW( J) ) 

TF(  1)=F(CL) 

DO  208  K=J,J3 

CL=C(SNEW(K) ) 

SIGC=SO(CL) 

ATMP=AT(TNEW(K) ) 

BTMP=BT(TNEW(K) ) 

SMT=ST(TNEW(K) ) 

SIGT(K)=SGT(SMT) 

208  CONTINUE 
IF(MELT,EQol)GO  TO  212 
J2=J3+1 

209  DO  210  K=1,MV 

IF(SIGT( J2-K).GE.SIGT( J3-K) ) GO  TO  210 

L=J2-K 

CALL  OVRTRN  (J,L) 

GO  TO  209 

210  CONTINUE 
CL=C(SNEW( J) ) 
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TF( I)=F(CL) 

IF(TNEW(J),LE.TF( I) ) GO    TO    212 

IF(ZI ( I ).L6oERR(DTAU) ) GO    TO    212 

MELT=1 

DO    211    K=1,MV 

IF(SNEW( J+K).EQoSNEW(J) ) GO    TO    211 

JSA=K 

GO    TO    202 

211  CONTINUE 
JSA=MV1 
GO    TO    202 

212  OK  I)=HIF(  I  )/(  RHOI*HS) 
CHIFU  )  =  CH1F(  I)+HIF(  I)*DTAU 
HF(  I)=THF*ZKW 

CHF(  I)=CHF(  I)  +  HF(I )*DTAU 

AH( I)=AHF( l)-HF{ I) 

CAH( I)=CAH(  I  )+DTAU*AH( I ) 

HC( I)=HF< I)-HIF( I ) 

CHC( I)=CHF( I)-CHIF( I) 

213  CONTINUE 
I1=MV1+1 

DO  301  I=I1,MVH1 
T(I  )=TNEW( I ) 
S(  I)  =  SNEW(  I ) 
301  CONTINUE 
J=MV1+1 

DO  403  I=2,MH1 
RHO=OE(SIGT( J) ) 
CP  =  SH(S< J)  ) 
VEC=VTK*RHO*CP 

IF(ZI( I ).LE0ERR(DTAU) ) GO  TO  401 
TI=TF(I ) 
CI=CO+BTC*SI/TI 
GRADI=DTEMP/ZI ( I ) 
QI=CI*GRADI 
TGRADS=QI/VEC 

TGRADR=0.35*QR/VEC 
TGOW=QS/VEC 

IF(TGRADSoGT.TGOW)GO  TO  401 
GO  TO  402 

401  TGRADS=QS/VEC 
TGRADL=QL*EXP(-0.5*ZI( I) )/VEC 
TGRA0R=(0.  35*QR+0a65*QR*EXP(-C 

402  STGD( I  )=TGRADS+TGRADL+TGRADR 
HUF( I )=STGD( I )*VEC 
PUF( I )=HUF( I )*DTAU 
J=J+MV1 

403  CONTINUE 
IF(AM0D(TAU,36'JJo  ) 
GO  TO  101 

501  Jl=l 
TIME=TAU/3600a 
J2=MVH1-MV 

WRITE(6,200)TIME,IRUN 
WRITE ( 6, 25C ) (TF(K) ,K= It MH1 ) 
DO  502  1=1, ^Vl 
WRITF.<  6,300)(T(  J)  , 
WRITE(6,35,  ) (S( J), 
WRITE (6, 400) (SIGT( 
J1=J1+1 
J2=J2+1 

502  CONTINUE 
WRITE (6,450) 
DO  50  3  I=2,«H1 
WRITE(6,5CMHIF(I) 

1HF(I)  ,CHF(  I),AH( I) 

503  CONTINUE 
IF(TIME.GEo 
GO    TO    101 

100    F0RMATU2) 
150    FORMAT ( 11F7 


5*ZI(I) ) )/VEC 


LT. l.E-4)    GO    TO    501 


J=J1 ,J2,MV1) 
J=J1, J2,MV1) 
J ), J=J1, J2,MV1) 


,CHIF(I) ,DI (I ),ZI(I) ,HC(  I),CHC(I), 
,CAH(I) ,AHF( I) ,CAHF( I) 


48, )G0    TO    1 


2) 
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200  FORMAT  (  US  'PROFILE  ( X-Z  PLANE)  OF  LEAD  AT  TIME=f,F5.1 

1  ,  10X, 'RUN  NUMBER' ,13) 
2  50  FORMAT (  ■ 0«  , ' TF '  , 3X , 21 ( 1 X  ,F5. 2 )  ) 
300  FORMAT ( «C» ,'T' ,4X,21( 1X.F5.2) ) 
350  FORMAT (•  •  ,  • S ■  ,4X, 2 1 ( i X , F5. 2 )  Y 
400  FORMAT (•  ' , • ST ' , 3X , 21 ( 1 X, F 5. 2 ) ) 
450  FORMAT (//, ' G ' , 3X , » H I F • , 7 X, » CHI F ■ , 6X, • DI « , SX , « Z I ■ , 8X , 

l'HC  ,8X,  'CMC  ,  7X,«HF'  , 8X,*CHF«  ,7X,»AH'  ,8X,  'CAH'  ,7X, 

2«AHF« ,7X, 'CAHF', /) 
500  FORMAT (•  « , 1 PE9. 2, 1 1 ( 1 X , E9.2 ) ) 
10C0  STOP 
END 

SUBROUTINE  CVRTRN  (I,K) 

CCMMON  TNEW (232 ) , SNElH  232) ,SIGT( 232) 

C(A)  =  (  A-0.030)/1.8(  5(x 

S0(  A  )=-C .069+1.  47Q8*A-0«0C'l  570*  A**2+0o0GC 0398*  A**3 

AT ( A )  =  A* ( 4  o 786 l-r ■ . 098 1 8  5*A  +  0e 00 1 Q843*A** 2 ) *C . 001 

BT<A)=A*(18o03r-.?o8164*A+0.01667*A*#2)*Cc0CC0Gl 

ST(A)=-( A-3«98)**2/503.570*( A+283.  )  /(A+67.26) 

SGT(A)=A+(SIGO+G.1324)*(1.0-ATMP+(BTMP*( STGO-G. 1324) ) ) 

K1=K-1 

J$=l 

L=0 

DO  1  M=I ,K1 

L  =  L  +  1 

IF(SIGT(K-L)oGE.SIGT(K))GO  TO  1 

L1=1+(K-L) 

J$=2 

GO  TO  2 


1 

CONTINUE 

2 

CTEMP=C.O 

CSALT=O.C 

IFUS.EQ.l  )L1=I 

DO  3  N=  1. 1 ,  K 

CTEMP=CTEMP+TNEW(N5 

CSALT=CSALT  +  SNEVv'(N) 

CONTINUE 

NN=1+(K-L1 ) 

ATEMP=CTcMP/FLOAT( NN) 

ASALT=CSALT/FLOAT(NN) 

CL=G( ASALT) 

SIGO=SO(CL  ) 

ATMP=AT(ATFMP) 

BTMP=BT(ATEMP) 

SMT=ST(ATEMP) 

ASIGT=SGT(SMT) 

DO    4    M=L1,K 

TNEW(M)=ATEMP 

SNEW(M)=ASALT 

SIGT(M)=ASIGT 

CONTINUE 

RETURN 

END 
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